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Abstract: We present computer-assisted methods for analyzing stochastic models of gene reg- 
ulatory networks. The main idea that underlies this equation- free analysis is the design and 
execution of appropriately-initialized short bursts of stochastic simulations; the results of these 
are processed to estimate coarse-grained quantities of interest, such as mesoscopic transport 
coefficients. In particular, using a simple model of a genetic toggle switch, we illustrate the 
computation of an effective free energy $ and of a state-dependent effective diffusion coeffi- 
cient D that characterize an unavailable effective Fokker-Planck equation. Additionally we 
illustrate the linking of equation-free techniques with continuation methods for performing a 
form of stochastic "bifurcation analysis" ; estimation of mean switching times in the case of a 
bistable switch is also implemented in this equation-free context. The accuracy of our methods 
is tested by direct comparison with long-time stochastic simulations. This type of equation-free 
analysis appears to be a promising approach to computing features of the long-time, coarse- 
grained behavior of certain classes of complex stochastic models of gene regulatory networks, 
circumventing the need for long Monte Carlo simulations. 

1 Introduction 

Various ways to model gene-regulatory networks exist, ranging from logical (boolean), 
to stochastic (Monte Carlo methods) or deterministic (ordinary differential equations) 
models (for recent reviews see [UHl GUI ) • Each modeling approach has its advantages 
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and disadvantages. One advantage of stochastic modeling is that it takes into account 
fluctuations due to the inherently random nature of biochemical reactions. This intrinsic 
noise gives rise to significant effects when either the molecular abundances of protein 
or mRNA molecules are small or the kinetics of the transitions between the chemical 
states of the promoter are slow [2*31 |2~T| . 

The established approach for stochastic modeling of spatially homogeneous chemical 
systems was introduced by Gillespie [U"" • The Gillespie Stochastic Simulation Algorithm 
(SSA) is based on repeatedly answering two questions: when does the next chemical 
reaction occur and what kind of reaction is it? Gillespie [U"] derived a simple way to 
answer these two questions that reduces the problem to a continuous-time discrete space 
Markov process. 

The SSA generates exact sample paths of the stochastic process and, for sufficiently 
large networks, it is computationally more efficient than solving the chemical master 
equation. However, the large size of naturally occurring gene regulatory networks makes 
even the SSA computationally intensive and practically impossible to use for comput- 
ing the long-time behavior of the network. Consequently, an important restriction of 
stochastic computations for many networks of interest is that we can efficiently run 
stochastic Gillespie-based simulators for short times only. It is therefore natural to look 
for computational methods that use only short time simulations (and as few of these 
as necessary) to compute the required information for the system. Such a computer- 
assisted approach is presented in this paper. 

Model reduction often provides a natural path to efficient simulation of a compli- 
cated model. As in other branches of physical modeling, separation of time scales can 
lead to successful model reduction in gene regulatory network modeling. Separation of 
time scales is frequently present in this context because synthesis and degradation of 
new proteins and transcripts usually occurs on a slower time scale than processes that 
change the chemical state of proteins (e.g., mult imerizat ion, protein/DNA interactions, 
phosphorylation). Theoretical methods for stochastic model reduction that take advan- 
tage of separation of time scales are being developed (e.g. "2H1 13 HZ1 EH])- Analytical 
reduction techniques assume that fast variables are in quasi-steady state with respect to 
the remaining slow variables. If the quasi-steady state distributions conditioned on the 
slow variables can be determined, then they can be used to eliminate the fast variables. 

Our approach is also based on (and takes advantage of) the separation of time 
scales and the approximation (computationally, on the fly) of quasi-steady marginal 
distributions (conditioned on the slow variables). The main feature of our approach, as 
will become apparent through its description and illustration, is that we do not "first 
reduce and then simulate the reduced model" ; our methods come in the form of wrappers 
around a black box dynamic simulator, and could equally well be applied to the most 
detailed stochastic version of the network model or to its best explicit reduction already 
available. In our approach, results about the long-term dynamic behavior of a stochastic 
simulator do not come from long-term simulation; they come from the design, execution 
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and processing of the results of "intelligently designed" short bursts of direct dynamic 
simulation. 

We believe it is useful to draw here an analogy with the study of nonlinear dynamics 
in systems of ODEs. Long-term information in the form of detailed bifurcation diagrams 
can be obtained from long dynamic integration; yet the same information is much more 
systematically and economically obtained through different algorithms using the same 
model: bifurcation, stability and continuation methods. It is this alternative to direct, 
long-term stochastic simulation (whether with the full detailed network model or with 
any good analytical reduction of it) that our approach makes available to the mod- 
eler. Ours is a "design of computational experiments" approach; it is guided by model 
reduction, but a reduced model is never explicitly obtained. 

The remainder of the paper is organized as follows. In Section El we introduce the 
genetic toggle switch as a simple model to illustrate our methods, and we specify the 
main questions that one would like to answer with these techniques. In Section El we 
present the general mathematical framework and main ideas of equation-free analysis 
[23 El U21 EH EH] • In Section HI we present an analysis of a deterministic model of 
the genetic toggle switch to provide insight into this system. We also introduce several 
stochastic models of increasing complexity that are used to illustrate equation-free anal- 
ysis. In Section El we compute the effective free energies and the associated stationary 
distributions for the stochastic models described in Section |3J Equation-free bifurcation 
analysis is then presented, and, in bistable cases, the mean first passage times for the 
system to switch between apparent stable fixed points are computed. We end with a 
discussion of the equation-free approach, its strengths, weaknesses, relations to other 
current methods for the acceleration of SSA-type simulations (e.g. [2B1 El UZ1 EH E]) 
and its possible extensions in Section El In particular, we will discuss the applicability 
of our methods to more complicated gene-regulatory networks. 

2 Model Description 

Our illustrative example is a two gene network in which each protein represses the 
transcription of the other gene (mutual repression). This type of system has been 
engineered in E. coli and is often referred to as a genetic toggle switch (HI EE]- The 
advantage of this simple system is that it allows us to test the accuracy of equation-free 
methods by direct comparisons with results from long-time stochastic simulations. In 
Section El we discuss the applicability of our methods to more complex problems where 
long direct stochastic simulation is impossible and the accuracy must be checked by on 
line a posteriori error estimates. 

A simple version of the genetic toggle switch is schematically drawn in Figure ^ 
The system contains two proteins Pi and Pi- The production of Pi (P2) depends on the 
chemical state of the upstream operator 0\ (O2). If 0\ is empty then Pi is produced at 
the rate 71 and if 0\ is occupied by a dimer of P2, then protein Pi is produced at a rate 
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Figure 1: A schematic diagram of the genetic toggle switch. 



e\ < 7i. Similarly, if 2 is empty then P 2 is produced at the rate 7 2 and if 2 is occupied 
by a dimer of Pi, then protein P 2 is produced at a rate e 2 < 7 2 . Note that for simplicity, 
transcription and translation are described by a single rate constant. The biochemical 
reactions and rate constants that correspond to the processes shown in Figure Q are 
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where the overbars denote complexes. Equation (|2.1|) describes production and degra- 
dation of protein Pi, equation (|2.2jl describes production and degradation of protein 
P 2 , equations (|2.3|) and (|2.4|) are dimerization reactions and equations (|2.5|) and (|2.6|) 
represent the binding and dissociation of the dimer and DNA. 

Single cell fluorescence measurements can be used to measure intercellular variability 
in protein expression levels. Therefore it is important to have efficient methods for 
computing the steady-state distribution of protein abundances from stochastic models 
similar to the one defined by (j2.1j) - (j2.fi j) . For moderately complex systems using 
long-time Monte Carlo simulations quickly becomes computationally prohibitive. We 
will illustrate how equation-free analysis can overcome this difficulty by accelerating the 
exploration of certain features of the long-term dynamics of the stochastic simulation. 
For certain values of the model parameters, the genetic toggle switch is bistable. If the 
system is described in terms of ordinary differential equations (ODEs) for the protein 
concentrations, then standard bifurcation analysis (numerical continuation methods) 
can be applied to determine the regions of parameter space in which bistability occurs. 
Using the model described by (j2.1|) - ()2.6|) as an example, we show how to extend these 
techniques to stochastic models. An important quantity that characterizes the dynamics 
of bistable stochastic systems is the average time for spontaneous transitions between 
stable steady states to occur. We will illustrate how this mean first passage time can be 
computed by using only short-time simulations. 

3 Equation-Free Analysis: Mathematical Framework 

Let us suppose that we have a well-stirred mixture of chemically reacting species; our 
main assumption is that the evolution of the system can be described in terms of a single, 
slowly evolving random variable Q (the approach carries through for the case of a small 
number of slow variables, but in this paper we will focus on the single slow variable 
case). Q might be the concentration of one of the chemical species or some function of 
the concentrations. Let R denote a vector of the other (fast, "slaved" system variables). 
Our assumption implies that the evolution of the system can be approximately described 
by the time-dependent probability density function f(q,t) for the slow variable Q that 
evolves according to following effective Fokker-Planck equation [33J: 

% = 5(|W(«.')1-W('.4 (") 

If the effective drift V(q) and diffusion coefficient D(q) could be explicitly written down 
as function of q, then (|3.1|) can be used to compute interesting properties of the system 
(e. g., the steady state distribution). Note that in addition to the assumption of a 
single slow variable, the validity of equation (j3.1|) requires sufficiently large molecular 
abundances and sufficiently fast chemical kinetics for transitions in the chemical state 
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of the operator jniEBj- Assuming that (j3.1j) provides a good approximation, we make 
use of the following formulas for the drift and diffusion coefficient fHJ 1201 ESI EZ| 



= <Q( t + At)- q m) = q > 

W At^O At v ; 

= 1 <[Q(« + At )-,]»|QM=«> 
w 2 At^o At v ; 

As described below, estimates of these two quantities can be found by using short-time 
bursts of appropriately initialized stochastic simulations. The steady solution of ()3.1|) is 
proportional to exp[— /3$(g)], where the effective free energy $(g) is defined as 

$(g) _ f q V(q' 
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/3$(g) = - / — )^-&q' + \nD(q) + constant. (3.4) 
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Consequently, computing the effective free energy and the steady state probability distri- 
bution also can be accomplished without the need for long-time stochastic simulations. 
A procedure for computationally estimating V(q) and D(q) is as follows: 

(A) Given Q = q, approximate the conditional density P(r\Q = q) for the fast 
variables R. Details of this preparatory step are given below. 

(B) Use P(r\Q = q) from the step (A) to determine appropriate initial conditions 
for the short simulations and run multiple realizations for time At. Use the results 
of these simulations and the definitions ()3.2j) and ()3.3)1 to estimate the average 
velocity V(q) and effective diffusion coefficient D(q). 

(C) Repeat steps (A) and (B) for sufficiently many values of Q and then compute 
$(g) using formula ()3.4|) and numerical quadrature. 

A very important feature of this algorithm is that it is trivially parallelizable (different 
realizations of short simulations starting at "the same q" as well a simulation realizations 
starting at different q values can be run independently, on multiple processors). 

In order to use the algorithm (A) - (C), we have to specify how the step (A) is 
performed. There are several computational options to approximate the conditional 
density P{r\Q = q). The simplest approximation is to estimate (through numerical 
experiments) the conditional mean <R|Q = q> and approximate P(r\Q = q) as a Dirac 
delta function 5{r— <R|Q — q>). Then the step (A) reads as follows: 

(Al) Given Q = q, pick an initial guess for the conditional mean of R. Denote 
the initial guess as <R(0)>. Run multiple realizations for a short time 5t and 
compute <R(5t)> . This procedure defines the mapping <R(0) >^<R(5t) > . 
Find the steady state of this mapping using standard numerical methods. The 
steady state is the required conditional average <H\Q = q> . Initialize R(0) as 
<R|Q = q> in all realizations in part (B) of the algorithm. 
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Another option is to approximate P{r\Q = q) as a distribution characterized by a few 
parameters, e.g. as a Gaussian distribution with mean ^ and variance cr. This can be 
done as follows: 

(A2) Given Q = q, pick initial guesses for the mean /x(0) and variance <x(0) 
of the conditional distribution function P(r\Q = q). Use this distribution to 
generate many realizations of R(0). Using these realizations as initial conditions, 
run stochastic simulations for a short time St and compute R(<5t). Computing mean 
and variance of R(<5t), we obtain the mapping [n(0), cr(0)] — > [/-i(5t), er(5t)}. Next 
use standard numerical methods to find the steady state [ft, cr] of this mapping 
and approximate P(r\Q = q) as a Gaussian distribution with mean and variance 
cr. 

The conditional density P{r\Q = q) can be also approximated by other basis func- 
tions. It is straightforward to generalize (Al) or (A2) to such a case. The better the 
approximation of P(r\Q = q) we have, the shorter the time step, 5t, required in the 
step (B) to achieve the same accuracy. So, a better approximation of P{r\Q = q) in 
step (A) decreases the computational intensity of step (B). On the other hand, step 
(A) is more computationally intensive if we want to obtain a better approximation of 
P(r\Q = q). One possibility for generating a better approximation of P(r\Q = q) is to 
use a "run-and-reset" procedure as was done in [Hj. This is accomplished as follows. 

(A3) Given Q = q, initialize the other variables R = R(0) of the system. Run 
stochastic simulations for the short time 5t. Then reset the value of Q(St) to its 
original value q keeping R unchanged. Repeat this procedure for many time steps 
and compute the conditional density P(r\Q = q) as a histogram of the recorded 
values of R. 

The approach (A3) attempts to compute the P(r\Q = q) effectively by successive sub- 
stitution, without resorting to numerical algorithms of the Newton-Raphson type for 
locating fixed points of mappings; we will return to this latter issue in the Discussion 
section. In our illustrative computations in Sectional we use step (A) in the form (Al) 
or (A3) for the simple stochastic models described below. Both give good results for our 
illustrative example. Since (Al) works for sufficiently long times 5t, there is no need to 
use (A2) or higher order approximations. For some stochastic simulations of our model 
problem, we also use slightly modified versions of the methods (Al) or (A3) as will be 
described in Section 

3.1 Bifurcations 

In deterministic problems, we often summarize the parametric dependence of the long- 
term dynamics in terms of bifurcation diagrams; for example, we may plot the steady 
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states of a deterministic set of ODEs as a function of a distinguished bifurcation param- 
eter. Several excellent continuation methods have been developed, implemented and 
made available for this purpose over the years, such as AUTO [Hl[7j. 

Here we illustrate how these methods can be extended to stochastic models [2U 
H2J EH I2HJ 123 • We assume, as above, that we have a stochastic problem that can be 
effectively described by a single variable Q. Let 7 be the bifurcation parameter. The 
first two steps in the algorithm are as follows: 

(A) Given Q = q and the value of the bifurcation parameter 7, compute the 
conditional density P{r\Q = q) using step (A) of the previous algorithm. 

(B) Using P(r\Q = q) from step (A) to determine the initial conditions, run 
multiple stochastic simulations for a short time At and compute the conditional 
average <Q(At)\Q(0) = q>. 

Steps A and B define the mapping (Q(0), 7) — ><Q(At)>. We denote this mapping as F, 
i.e. F(Q,j) —<Q(At)> . Our goal is to track the fixed points of F (i.e. F(Q,j) = Q) 
as the bifurcation parameter 7 is varied. To do this, we first use a Newton-Raphson 
algorithm to find two steady states (Qi,7i) and ((52,72) which are sufficiently close 
to each other (note that one can estimate the derivative of F{Q 1 r y) numerically by 
evaluating F{Q ) r f) at different points). Then, in a parameter continuation context, we 
choose a small parameter 5 (which can be modified adaptively during the computation) 
and find the next steady state using continuation. That is, the next values of Q and 7 
are found by solving the following system of equations 



To find the solution of (J3.5|) . we estimate the Jacobian numerically by evaluating F(Q, 7) 
at several points and then use Newton-Raphson algorithm. When the number of vari- 
ables starts becoming large, matrix-free methods of iterative numerical linear algebra 
(such as Broyden, or Newton-Krylov GMRES [221 ) can be used to solve for the fixed 
point, as opposed to full numerical Jacobian estimation. The fixed points computed 
this way provide, under certain conditions, good estimates of the critical points (min- 
ima, saddles) of the effective potential as a function of a model parameter 7; this 
issue is discussed extensively in [2E1 123 El E3 , and we will return to it again in the 
Discussion section. 

3.2 First Passage Time 

Suppose that we have a bistable stochastic system. That is, the effective free energy $(g) 
has two local minima [14] - see Figure 121 Then an important quantity characterizing 
the long-time system dynamics is the mean time for spontaneous transitions to occur 
between the stable steady states. Let q m < qu denote the two stable steady states and 
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Figure 2: Potential 0/ £/te bistable system. 



let be the unstable state (i.e. local maximum of Then we define the first 

passage time for transitions from q m to qu as 2r e where r e is the average time for the 
system to reach the unstable steady state q u for the first time given that it starts at q m . 
The factor of two occurs because once the system reaches the unstable steady state, half 
the time it returns the original stable steady state q m and the other half of the time it 
transitions to qM ■ 

Algorithm (A) - (C) gives a procedure to estimate the effective potential by 
running short simulations only. Once we have the effective potential, we can compute 
To as follows 
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Equation (J3.6)) can be further simplified if the height of the potential barrier [&(q u ) — 
$(g m )] is large compared to the noise strength. In this case, the limit q in the second 
integral can be replaced by q s , allowing the two integrals to be evaluated separately 
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The main contribution of the first integral stems from the region around q s , and the main 
contribution from the second integral stems from the region around q m . Consequently, 
we expand <3>(g) according to 
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for the first and the second integral, respectively [33]. When these expansions are used 
in equation ()3.2jl . the following result is obtained 

ATik B T 

~ [D{q u ) + D(q m )W&'(q m )\V'(q u f\ 

which is the generalization of Kramers' formula to the case of a state dependent diffusion 
coefficient [331 EH|- Formulas (|3.6|) and (|3.7|) are both used in Section 15^1 to estimate r e . 
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4 Analysis of the Model Problem 

In this section, we study the behavior of the model given by equations (j2.1|) - (|2.6|) . 
To provide insight into the problem, we start by analyzing the deterministic system. In 
Sections 14.21 and 14 .3( we introduce two stochastic models that are simplified versions 
of the model defined by (|2.1|) - (|2.6|) . We use these models because of the relative 
ease in performing long-time stochastic simulations with them; this allows the results 
from the equation-free analysis to be validated by direct comparisons with Monte Carlo 
simulations. We will also verify that the equation-free methods can be applied to the 
full model. As discussed below, for this case the long-time Monte Carlo simulations 
become computationally very expensive. 



4.1 The Deterministic Model 

To simplify the deterministic analysis, we make the assumption that equations ()2.3j) - 
(I2.fjjl are at quasi-equilibrium and derive deterministic rate equations for the protein 
concentrations. Let x\ and x 2 denote the average monomer concentrations of P\ and 
P 2 , respectively, and let d\ and d 2 denote the respective dimer concentrations. Also, let 
o\ and 02 denote the probabilities that the operators 0\ and 2 are not occupied. For 
the dimerization process the assumption of quasi-equilibrium implies 

d\ = - — x 1} and a 2 = 7 — x 2 . (4.1) 

K-i K_2 

Similarly, the quasi-equilibrium assumption for the operators implies that 

°i = I k ~°) , , and o 2 = k ~° 2 . (4.2) 

fc_oi + k ol d 2 k_ o2 + k o2 di 

The total concentration of P\ is given by y\ — x\ + 2d\. The total concentration y\ 
evolves according to the following ordinary differential equation 

^ = 7i0i + - ox) - 5ix\. (4.3) 
at 
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where 5 is the degradation rate of the monomers, and it has been assumed that dimers 
are protected from degradation. Substituting y± = x± + 2d\ — X\ + 2-^xf into (|4.3|) . we 
obtain 



. ki \ dxi 

Finally, using (|4.1|) - (|4.2|) produces 
dxi 1 
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where the parameters K\ and u)\ are defined as follows 
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Using similar reasoning an analogous equation for x 2 can be derived. 

For simplicity, we will present the symmetric case in which the rate constants for 
processes involving Pi are identical to those involving P 2 . That is, we assume k = 
K i — K 2, 7 = 7i = 72, ^ = o»i = u 2 , and 5 = 5i = 5 2 . Moreover, we assume that 
the production rate is zero if an operator is occupied, i.e. E\ = e 2 = 0. Making these 
assumptions, (|4.4jl simplifies to 



dx\ 
~dF 



1 + KX\ 
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(4.6) 



and the equation for x 2 is obtained by alternating the subscripts in the above equation. 
Hence, the problem has been reduced to a system of two equations with four parameters. 
Note that the value of k does not influence the steady-state behavior of the system. In 
this paper, we fix the values of 5 and to to be 0.00075 and 2 x 10~ 6 , respectively. 

The steady states values of function of 7 are shown in Figure El In this 

figure, solid lines denote stable steady states and dashed lines denote unstable steady 
states. For 7 < 1.06 there is a single steady state. At 7 = 1.06 a pitchfork bifurcation 
occurs, and for 7 > 1.06, there exist three steady states. The steady state with 37 = x 2 
is unstable and the other two steady states are stable. 

Due to separation of time scales, the long-term dynamics of this problem lie on a 
lower- dimensional (here one-dimensional) slow manifold; this suggests that one may be 
able to construct an effective one- dimensional dynamical system describing the long- 
term evolution of the model on (near) this slow manifold. In constructing such a re- 
duced model, an important question even in the simple deterministic case is the choice 
of the right observable - the variable in terms of which the long-term dynamics will be 
expressed. An extensive discussion of the choice of such a "right observable" for the 
deterministic case can be found, for example, in [11] ; as discussed there, even if we do 
not know the exact slow variables, any set of variables that parametrizes the slow man- 
ifold can be practically used to reduce the system in an equation-free context. For the 
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Figure 3: The dependence of the steady state values of xi on 7. The solid lines denote 
stable fixed points and the dashed line corresponds to unstable fixed points. In this figure 
and throughout the paper 5 = 0.00075 and uj = 2 x 10~ 6 . 

stochastic case, a good early illustration and discussion of manifold parametrization can 
be found in [36] . Choosing the right observable is an important issue in the implemen- 
tation of equation-free computations, and the subject of intense current research which 
we will briefly comment on in Sectional 

In this paper, and for this example, our equation-free analysis assumes that the 
problem can be described in terms of a single variable. Consequently, it becomes im- 
portant to select a good observable that further simplifies the two-dimensional problem 
to one dimension. A tempting (and obvious) choice for the one-dimensional observable 
is the molecular abundance of P\ (or P2). We demonstrate below that using Pi in the 
equation-free analysis produces good results. However, we also make use of the symmet- 
ric variable defined as the difference in the protein abundances Q = Pi — P%. In terms 
of the rate equations the symmetric variable is s = X\ — x<i- The bifurcation diagram in 
terms of s is shown in Figure UJ The symmetry of the diagram suggests that Q might 
be a more natural observable than Pi (which also produces good results, as we will see 
below). 

4.2 Stochastic Model I 

To start our investigations in the equation-free framework, we constructed a very simple 
stochastic model of the system. We use this simple model to benchmark equation-free 
computations, since the results can be tested against Monte-Carlo simulations easily. 
Results for the full system are also presented below. The simple stochastic model consists 
only of reactions for the synthesis and degradation of proteins Pi and P2, but the 
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Figure 4: The dependence of steady state values of the symmetric variable s = X\ — %2 
on 7. 



following effective rate constants are used 

1 7 

l + K-Pl l + ^Pj 

^ Pi (4.7) 

1+rePj 



1 7 
1 + kP 2 1+wP-j 



S v ' 

1 + kP 2 

The above reactions are consistent with the deterministic model, but in general do not 
preserve the noise structure of the full stochastic model. 

To simulate the mechanism contained in model (|4.7jl - ()4.8jl . we use the standard 
Gillespie SSA [IB]. The results for different values of the parameter 7 are plotted in 
Figure El For each 7 we plot the time evolution Pi (left panel) and Q = Pi — Pi (right 
panel). We see that for small 7, the solution fluctuates around the stable deterministic 
steady state with relatively small noise amplitude. When 7 = 1.06, the noise amplitude 
has increased substantially, which is typical of stochastic systems near a "bifurcation" ; 
the word bifurcation is put here in quotes to denote that (in contrast to the deterministic 
case) there is no isolated parameter value marking the onset of bistability - no clear 
bifurcation point exists for the stochastic dynamics. Yet one can still claim that a clear 
bifurcation point exists for the critical points of the potential $(5,7) in the stochastic 
model; furthermore, depending on the time horizon of our observation of a stochastic 
simulation, one may still appear to see an apparent bifurcation point for its averaged 
statistics (see the discussion in [TBIE])- If 7 is increased further, then Q = is no longer 
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Figure 5: Stochastic Model I. Plots of P± and Q = P\ — P 2 (is a function of time for 
different values 0/7. The parameter values used to produce these figures are 5 = 0.00075, 
u = 2 x KT 6 , and k = 2 x 10~ 4 . 
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a stable steady state and the system clearly shows bistability. All plots are computed for 
the same time interval [0, 15 x 10 6 ]. For 7 = 1.25 the steady states are sufficiently stable 
so that no transitions occurred in this time interval (data not shown). Therefore, for this 
case, determining the steady state probability distribution from long-term Monte-Carlo 
simulations would be very time consuming. 



4.3 Stochastic Model II 

Stochastic Model I considers only two variables P\ and P^- Here, we introduce a stochas- 
tic model that also takes into account the biochemical states of the operators, while 
maintaining the assumption that the dimerization reactions (12. 3j) - ()2.4j) are at equilib- 
rium. That is, we consider the four variables Pi, P2, 0\ and O2. The model is defined 
in terms of the following reaction steps: 



l + K-Pl 



p i (4.9) 



— 

s 

1+K-P2 
K 

u o 1 = o" 

K 

"o 2 = o" r~ ' 

and contains an extra parameter, K = fc_ ol = k_ o2 . Note that u O\ = 0" means that the 
operator 0\ has a dimer of P2 bound to it and therefore is "off" and "Oi = 1" means 
that the operator 0\ is empty and therefore "on". The same is true for O2. This implies 
that the random variables 0\ and O2 are binary, whereas the variables Pi and P2 can 
take on any non-negative integer value. Stochastic Model I is recovered from Stochastic 
Model II in the limit K — > 00. We thus expect the models to produce similar results for 
large values of K. 

Again, we use the standard Gillespie SSA [THj to simulate model (14. 9 J) - (|4.12|) . 
The results for different values of K for 7 = 1.14 are plotted in Figure El Comparing 
Figure El and corresponding panel from Figure El we can confirm that Stochastic Model 
II produces the same behavior as Stochastic Model I for large K. However, in general, 
different values of K can change the bifurcation structure of the system and affect the 
first passage times between the two stable steady states of the bistable system [2~3] . 

Because Stochastic Models I and II do not explicitly take into account dimerization, 
which in general is a fast process, they run much more efficiently than the full model 
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Figure 6: Stochastic Model II. Plots of P\ and Q = P± — P2 as a function of time for 
different values of K and 7 = 1.14. The other model parameter values are the same as 
in Figured 
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given by (j2.1j) - (|2.fij) . However, they do not in general preserve the noise structure of the 
full system. In the next section we use all three models to highlight the computational 
features (and potential benefits) of equation-free analysis. 

5 Results of Equation-Free Analysis 

In our approach, we want to study the stochastic models presented above using only 
short bursts of appropriately initialized stochastic simulations; the goal is to design these 
bursts, and process their results so as to determine long-time properties of the system 
(e.g., steady-state distributions, bifurcations, mean first passage times) efficiently We 
use (and compare) the different algorithms discussed in Sectional 

5.1 The effective potential and steady state distribution 

In this section, we use equation-free analysis to evaluate the effective potential (an 
"effective free energy" ) and the steady-state distribution for Stochastic Models I and II 
and the full system. We start with Stochastic Model I. First, we will consider the slow 
variable Q = P\ — Pi and the fast (slaved) variable R = Pi + P2- Initially the preparatory 
step (A) of the algorithm presented in Section El was done using the method outlined in 
(Al), i.e. we used the conditional mean <R\Q = q> to initialize the computations in 
the step (B). A good approximation to this average can be found using the deterministic 
equations, and this was the number used in our preliminary computations to initialize 
the simulations in the step (B). That is, for a given Q, we initialized all realizations 
in the step (B) with the same value of R. Then we chose At equal to 100 time steps 
of the Gillespie SSA. Note that this implies that the actual value of At varies for each 
realization and depends on the values of the rate constants. However, the computer 
(CPU) time is the same for all the results presented for this case. 

The equation-free results for the effective potential for different values of 7 are given 
in Figure [7| These results are in good agreement with the long-term stochastic sim- 
ulations presented in Section 14.21 The potential has a single minimum 7 < 1.06. As 
7 is increased the potential broadens implying the system becomes "noisier". When 
7 > 1.06, the potential shows two local minima and the system is bistable. 

Since we are using a very simple stochastic model, it is not computationally expensive 
to compute the steady-state distributions directly by long time simulations. We use the 
Gillespie SSA to generate 10 11 time steps of the stochastic process and recorded the value 
of Q at each time step. The resulting time series was binned to produce the steady- 
state distribution of the system. Figure |S] presents a comparison of the two computed 
steady state distributions. The results obtained by long-time simulations are shown as 
blue histograms and the steady state distributions computed from the effective potential 
C exp[— /3$(Q)] are given by the red lines. We see that equation-free analysis gives very 
good results. 
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Figure 7: The effective free energy $ for different values 0/7 computed by our procedure. 
The other model parameter values are the same as in Figured 

In Sectional we introduced three possible methods, (Al) - (A3), to perform the 
preparatory step (A) (typically called the "lifting" step in the equation- free framework). 
We have shown that approach (Al) produces good results for Stochastic Model 1. Since 
(Al) works, there is no need to improve the results by considering (A2). Instead, we 
discuss the (A3) approach. In this approach given Q = q we run the simulations for a 
short time 5t and record the value of R. Then we reset Q = q but leave R unchanged. 
We repeat this procedure many times and compute the conditional density P(r\Q = q) 
as a histogram of recorded values of R. In our simulations, we chose 5t equal to one SSA 
step. To compute the conditional density P{r\Q = q), we used 11 million SSA steps. 
First we let the system run for a million time steps to remove the transient in R, and 
then used remaining 10 million time steps to compute the conditional density. We used 
200 million Gillespie time steps in part (B) of the algorithm. Consequently, step (A3) 
did not significantly change the computational cost of the program. 

The graphs of P(r\Q = q) for 7 = 0.98 and 7 = 1.14 are given in Figure El The 
left panel in these figures shows P(r\Q = q) for five values of Q. The right panels 
show P{r\Q = q) as a function of r and q. Next, we can use the computed conditional 
density P(r\Q = q) to initialize R in the step (B). Doing this, produces results which 
are virtually identical to results from Figure |H] (graphs not shown). 

We now repeat the previous computations using the more complicated Stochastic 
Model II. The results are shown in Figures ITU1 and ITT1 In Figure ITUl we choose 7 = 1.14 
and compute the steady state distribution for Q for three values of K. The results 
are compared with direct simulations of Stochastic Model II and with each other. The 
results from Figure El can also be compared to the corresponding plot with 7 = 1.14 in 
Figure |Hl which can be viewed as the limit K — > 00. As can be seen, the results given by 
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Figure 8: Comparison of steady state distributions obtained from the effective free ener- 
gies shown in Figure^ with histograms (blue) obtained by long time simulations. 

Stochastic Model II for K = 10 are already in good agreement with the corresponding 
results obtained by Stochastic Model I. Figure ^2 shows similar results for 7 = 1.20. 
Again, we obtained accurate results using the equation-free method. 

Up to now we have used the symmetric variable Q = Pi — P 2 as our observable. 
However, we often do not have a priori knowledge of the slow variable, or, more generally, 
of a good observable to parametrize the long-time system dynamics. To investigate the 
sensitivity of our results to the choice of observable, we repeated the computations on 
Stochastic Model I using Pi instead of Q. To use Pi as our observable, we modify 
step (Al) so that we simply initialize P2 using P 2 = 1 1+ ^p^ ■ The numerical results for 
different values of 7 are given in Figure ED Again good agreement is seen between the 
equation-free method and the Monte-Carlo simulations. Because Pi has both a slow 
and a fast component, this result illustrates that equation-free methods can be used 
even when the slow variable is unknown. An extensive discussion of this point in a 
deterministic context can be found in jTTj: one does not necessarily need the correct 
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Figure 9: Conditional distribution P(r\Q = q) for Stochastic Model I. Pictures on the 
left show P(r\Q = q) for selected values of q. Pictures on the right show P(r\Q = q) as 
a function of r and q. 



slow variable - one needs an observable that parametrizes the slow manifold, a quantity 
in terms of which the slow manifold can be expressed as the graph of a function. 

Encouraged by the success of our computational framework for the simple Stochastic 
Models I and II considered above, we next investigated how well these methods would 
work on the full system described by equations (j2.1|) - (|2.6|) . We first performed long 
time Monte Carlo simulations using BioNetS [I]. A two dimensional histogram for the 
total protein numbers I\ = P\ + 1P\P\ and T 2 = P2 + 2P 2 P 2 is shown in Figure IT37 a) . 
This simulation took over 500 total CPU hours and is the result of 800 runs distributed 
over 18 CPUs (over a trillion Gillespie SSA steps). The black curve in Figure EHb) 
is the projection of the histogram onto the T\ axis. We next performed equation-free 
computations for the system. As our single observable, Q, we used the total protein 
number Ti, because this is a quantity that can be measured using single cell fluorescent 
techniques. We used a slightly modified version of step (A3) to compute the conditional 
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Figure 10: Comparison of steady state distributions from Stochastic Model II. The top 
left panel are results from the equation-free analysis for K = 0.1, K = 1 and K = 10 
and 7 = 1.14. The remaining three panels compare these results (red lines) to the steady 
state distribution computed from long-time Monte Carlo simulation (blue histograms). 
The other model parameter values are the same as in Figure^ 

density P(r\Ti = ti). For a given value of T 1; we set the rate constants for synthesis 
and degradation of this protein equal to zero. We then ran the simulations for a time 
of 1 x 10 5 to remove any transients. Next still keeping Ti fixed, 10000 samples of the 
other variables were collected at evenly space intervals over a time period of 2 x 10 5 
and used to generate the conditional density. A time step of At = 15 was used in step 
(B) of the algorithm. To compute the steady state distribution, polynomials were fit to 
the average velocity and effective diffusion coefficient computed from the equation-free 
analysis and then used to compute the effective free energy. The red curve shown in 
Figure fIBIf b) is the result of the equation-free analysis. It took less than an hour of 
CPU time. Very good agreement between the equation-free method and Monte Carlo 
simulation is seen. Our investigations into these methods revealed that whereas the 
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Figure 11: Stochastic Model II. Comparison of steady state distributions obtained by 
equation-free analysis (red line) with histograms (blue) obtained by long time stochastic 
simulations. In this figure 7 = 1.2 and K = 0.1 in the left panel and K = 1 in the right. 
The other model parameter values are the same as in Figure^ 



velocity, V(q), is robust to changes in At, the effective diffusion coefficient, D(q), is 
quite sensitive and needs to be treated with care. Also, because of the exponential in 
the integral for the effective potential, small changes in the average velocity or effective 
diffusion coefficient can have large effects on the steady state distribution. Therefore, 
it is important to average over sufficiently many realizations to ensure convergence of 
average velocity and effective diffusion coefficient. Better estimation techniques, such 
as those developed by Ai't-Sahalia using maximum likelihood [2] should be incorporated 
in the data processing step of the algorithms. Even with these caveats, the results 
presented in this section demonstrate the feasibility and high potential of equation-free 
methods for analyzing stochastic models of genetic networks. 



5.2 First Passage Time 

When 7 is sufficiently large the system is bistable. An important characterization of 
bistable systems is the average time for noise-induced transitions between the stable 
states. Here we make use of the definition of the first passage time from Section 1^21 For 
the results presented in this section we use Q = Pi — P2 as our observable and Stochastic 
Model I. The system is bistable for 7 > 1.06. Let the deterministic stable steady states 
of Pi be denoted as p m and pu with p m < pm- Because of the symmetry of our problem, 
p m and pm are also the stable steady states of P2. Let the random variable T e be defined 
as the first time when "Pi = P2" given the initial conditions Pi = p m and P2 = pu- In 
terms of Q, this means that T e denotes the time to reach Q = when the process starts 
with Q equal to the negative steady state q m = p m — pu- Let r e denote the average of 
T e . Then, direct Monte Carlo simulations can be used to compute the value of r e . The 
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Figure 12: Comparison of steady state distributions using the variable Pi as the observ- 
able for various values 0/7. The other model parameter values are the same as in Figure 
[3J Again, the red lines are the results of equation-free analysis and the blue histograms 
are obtained by the long-time stochastic simulations. 
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Figure 13: (a) The steady-state distribution for the total protein numbers computed from 
long-time Monte Carlo simulations of the full model (\2. 1\ - (jg. fij) . (b) The projection 
of the 2D distribution onto the T\ axis (black curve). The red curve is the result of the 
equation-free analysis. The parameter values used to compute these figures are 71 = 
72 = 1.14, 5 X = 5 2 = 7.5 x 10 4 ; e x = e 2 = 0, k x = k 2 = 5 x 1(T 4 , k_ x = k_ 2 = I, 
k i = k o2 = 0.004 and k_ ol = k_ Q2 = 0.1. These values are consistent with the parameter 
values K — 0.1, 5 = 7.5 x 10 4 and uj = 2 x 10~ 6 used in Stochastic Models I and II. 
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Pm 


Pm 


q m = Pm - Pm 


computed r e from simulations 


N 


1.14 


481.1 


1038.6 


-557.5 


7.0 x 10 5 ±6.7 x 10 3 


10 000 


1.20 


425.8 


1174.2 


-748.4 


1.6 x 10 7 ± 1.6 x 10 5 


10 000 


1.25 


392.4 


1274.3 


-881.9 


1.0 x 10 9 ±6.3 x 10 7 


250 



Table 1: The mean first passage time computed from long-time stochastic simulations, 
averaging over N transitions. The results are expressed in the form ([sample mean] ± 
[sample variance] /\f~N). 



results of such simulations for three different values of 7 are presented in the Tabled 
As expected, the computational time needed to compute the mean first passage time 
increases rapidly with 7. In Section 13 .2\ we introduced two formulas ()3.6|) and ()3.7|) 
to compute r e . Both formulas make use of the effective free energy computed by the 
equation-free algorithm. These potentials for 7 = 1.14, 7 = 1.20 and 7 = 1.25 are given 
in Figure [7| Consequently, we can compare the results obtained by the long simulations 
with the results found from formulas ()3.6|) and (j3.7J) for r e;p and r e;fc , respectively. The 
results are shown in Table EJ 

Not surprisingly, the results given by r e;p are better than results given by the Kramers' 
approximation r e -k- However both methods produce results that are within a factor of 
2 of the waiting times estimated from Monte Carlo simulations. As 7 becomes large 
the Monte Carlo simulations become computationally expensive. Therefore only 250 
realizations were used to estimate the mean first passage time, and we expect that the 
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r e from Table [T] 


r e . p given by ([3~T5]) 


r e . k given by ([3~T7]) 


1.14 


7.0 x 10 5 


6.1 x 10 5 


1.3 x 10 6 


1.20 


1.6 x 10 7 


1.4 x 10 7 


2.6 x 10 7 


1.25 


1.0 x 10 9 


6.7 x 10 8 


1.2 x 10 9 



Table 2: Comparison of the mean first passage time computed from equation-free analysis 
with long-time stochastic simulations. 

discrepancy between the Monte Carlo simulations and equation-free analysis for this case 
is due to finite sampling errors. Initializing the simulation at conditions that are rarely 
visited by the direct simulation itself constitutes a form of bias; this bias is designed to 
give faster computational estimates of the effective potential and - through this - of the 
first passage times. Clearly, this approach hinges on knowledge of a good observable, 
and in principle does not depend strongly on the value of the parameter 7; therefore, 
the larger the parameter 7 the higher the computational speedup in the first passage 
time estimation that will result. A quantitative study of this speedup is underway and 
will be reported elsewhere; it does not lie within the scope of this paper. We stress, 
however, that (as in molecular dynamics simulations) knowledge of a good observable 
(a good "reaction coordinate") is crucial for the success of the approach. 

Note that formula r e -k requires estimates of the second derivative of the potential at 
points q u and q m . To do this, we fit $(g) locally to a polynomial and used the derivatives 
of the polynomial at the required points; once again, maximum likelihood techniques 
(e.g. j2J) should be used for better results. The formula for r e;p , requires the evaluation 
of an indefinite integral. The integral was approximated by considering only a finite 
interval that neglected contributions from the region of sufficiently small q where the 
potential $ is very large. 

5.3 Bifurcations 

In this section, our goal is to run the simulations for short times only and compute a 
form of "stochastic bifurcation diagram" using continuation methods, as an extension 
of the deterministic bifurcation computations. We use Stochastic Model I and study 
the dependence of the "steady states" on 7; the "steady states" we report are the 
fixed points of the algorithm from Section I3~T1 with the conditional density P(R\Q = q) 
approximated by the Dirac delta function in (A) - (B), similar to the approach (Al) 
from Section |HJ Numerical results are given in Figure IT41 For comparison we also plot 
the steady states of the corresponding deterministic equation (compare with Figure EJ). 
The plot in Figure El was computed by initializing on different branches far from the 
bifurcation point and continuing from these different initializations (our simple arclength 
continuation algorithm did not include a "pitchfork detection" component). 

The accuracy of the numerical results depend on several factors: the estimation tech- 
nique for the Jacobian elements, the tolerance of the error for Newton-Raphson itera- 
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Figure 14: A plot of the steady states obtained by equation-free analysis fl,?.,5jl (blue 
circles). Also shown are the deterministic steady states from Figure^ (red line). 

tions, the number of realizations which are used to evaluate F, the time interval At and 
the steepness of the underlying potential $. As can be seen in Figure stochasticity 
along with all these numerical factors have slightly perturbed the pitchfork bifurcation; 
this could be exacerbated by our choice of (symmetric or asymmetric) observable. It is 
easy to follow any branch of steady states far from the bifurcation point. For obvious 
reasons this becomes more complicated when we are close to the "bifurcation point" at 
7 = 1.06. The main problem is that the potential becomes "flat" close to the bifurca- 
tion point — see Figure One way to improve the results is to adaptively change the 
number of realizations in (B) . That is, if the Newton- Raphson iterations of ()3.5)1 do not 
converge to a desired tolerance, then more realizations are added. Another approach is 
to estimate directly a local polynomial model of the underlying diffusion process from 
discrete SSA data using maximum likelihood tools, and then search for the bifurcations 
of the critical points of the effective potential. Indeed, one can plot the zeroes of the 
estimated drift, or - in the case of a state-dependent diffusion coefficient - one can cor- 
rect them to report the maxima of the steady state distribution [213 EZj; both of these 
are good candidate bifurcation diagrams for the stochastic case. When the potential is 
steep and the equilibrium is "less noisy" it is not necessary to use many realizations; the 
relation between computational effort (in terms of number of replicas, simulation time 
horizon and estimation method) and resulting accuracy is, again, a subject of current 
investigation beyond the scope of this paper. 

Finally, the results using Pi instead of Q as the observable are shown in Figure fTol 
In this case, the asymmetry of our observable, and the perturbation it causes on the 
initialization process, make the perturbation of the pitchfork bifurcation stronger. Of 
course, the results depend on the initialization procedure, our estimation technique, the 
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Figure 15: Plot of the steady states obtained by using P\ as the observable (blue circles). 
Again for comparison the deterministic case is shown as the red line. 

error tolerance, the number of realizations, the length of time step At as well as the 
type of continuation algorithm we are using (here we used a very simple one, without 
bifurcation detection, in order to demonstrate what is possible). Accurate bifurcation 
detection depends on accurate Jacobians and even higher derivatives; estimating these 
from dynamic (and noisy !) data is notoriously difficult. While conceptually we do 
have the tools to "hone in" the more accurate detection of bifurcation points, careful 
quantitative work is necessary to pin down the tradeoffs between computational effort, 
model estimation accuracy and bifurcation point estimation accuracy. 

6 Discussion 

In this paper we discussed and illustrated the use of certain equation-free numerical tech- 
niques that have the potential to accelerate the computer-assisted analysis of stochastic 
models of regulatory networks. There is a clear current need for accelerating such simula- 
tions: even for modestly complex regulatory networks, stochastic models rapidly become 
computationally expensive. Computational acceleration is usually based on model re- 
duction; theoretical methods for stochastic model reduction that take advantage of a 
separation of time scales are the focus of intense current research |2*H1 El El EH EI] • As 
we discussed in the introduction, many important gene regulatory networks do satisfy 
this assumption of a separation of time scales because synthesis and degradation of new 
proteins and transcripts usually occurs on a slower time scale than processes that change 
the chemical state of proteins. Analytical model reduction techniques assume that the 
fast variables are in quasi-steady state with respect to the slow variables, and use the 
quasi-steady state distributions conditioned on the slow variables to eliminate the fast 
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variables by averaging. These methods have been successfully applied to simple models, 
but the theory is not as well established as the deterministic counterpart. Having an 
explicit model lies often at the basis of such stochastic reduction methods. 

In the equation-free approach many of the same elements (separation of time scales, 
approximation of conditional quasi-steady distributions) also underpin computational 
efficiency; but the basic premise is that the model is available in the form of a "black 
box" simulation code. We do not try to first reduce, and then simulate the reduced 
surrogate; we try to design the smallest number of "intelligent" short computational 
experiments with the full stochastic model to find the quantities of interest, whether 
these are steady state probability distributions, their maxima, or transition rates in 
the bistable case. In that sense, the approaches we described here do not hinge upon 
the "inner", detailed simulator as being a Gillespie SSA one - the methods are equally 
applicable to any "inner" simulator, stochastic or deterministic, as long as the main 
assumption of a low-dimensional effective stochastic model is a good one for the long- 
term system dynamics. Indeed, if another reduction method can be used to produce 
a good approximate dynamic simulator, our algorithms can be "wrapped" around this 
surrogate simulator rather than the full model for further acceleration. 

Another important point has to do with the type of computation we are interested in 
- do we want to accelerate the direct simulation of the model, or do we want to accelerate 
the computation of certain features of its long-term dynamics (e.g. of the maxima of the 
steady state distribution)? These latter quantities can be also obtained from long-term 
direct simulation, but one of the points that we want to stress is that we can link direct 
simulation to different numerical algorithms (such as contraction mappings, and con- 
tinuation methods) to obtain these quantities, often faster than with direct simulation 
alone. In the same way that bifurcation diagrams for dynamical systems are usually 
not computed through direct ODE integration, but through bifurcation algorithms, the 
parametric dependence of the long-term dynamics of stochastic models does not have to 
be computed through long-time direct simulation only. This "alternative" acceleration, 
not through accelerating the direct simulation itself, but through linking it to different 
numerical algorithms, lies at the basis of the equation-free framework. 

Having said this, we briefly mention that equation-free methods for accelerating the 
direct simulation itself also exist. Coarse projective integration which uses short bursts 
of direct simulation to estimate time derivatives of evolving probability densities, and 
then passes them to standard numerical integration algorithms, has been successfully 
used in many contexts [2U El EE E| Coarse projective integration has a strong relation 
to the direct simulation acceleration methods in [15] ; it has not been discussed in this 
paper, because we chose to focus on very long-term features of the network dynamics; 
it might interest the reader that the method can be used to also integrate backward in 
time, and solve "effective boundary value problems" to find "coarse" limit cycles [32] . 

In the equation-free methods for analyzing stochastic models of gene regulation that 
we discussed in this paper, we have tried to circumvent the difficulties encountered by 
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direct simulation (in this case SSA) through the design of short bursts of appropriately 
initialized computational experiments with the full simulator. In a sense, we "resign 
ourselves" to the fact that the direct simulator is expensive; we ask what is the shortest 
amount of running of this expensive direct simulator in order to obtain the quantities 
we are interested in. The "design of experiment" protocols are templated on tradi- 
tional continuum numerical methods, like the fixed point and continuation algorithms 
to compute bifurcation points, or quadrature to estimate Kramers' formula. The only 
difference is that the quantities (residuals, actions of Jacobians, values of the integrand) 
that are required for numerical computation are not given by a closed formula, but 
rather through direct numerical simulation of the full model and estimation. We reit- 
erate once more that these techniques can be wrapped around the full direct simulator, 
or our best available reduction of it, without change. 

Knowing appropriate coarse-grained observables (the variables in terms of which the 
unavailable effective model would be written) is an important feature of the algorithms. 
Extensive experience with the problem, intuition, or analytical work may often suggest 
such observables; we did take advantage of such knowledge in this paper. We did 
already demonstrate an important point: more than one observables are capable of 
doing a good job as the parameterizing variables in an equation-free context; one does 
not need to know the exact slow variables. This issue is discussed extensively for the 
deterministic context in [TT]. It is, however, important to note that algorithms for the 
detection of low- dimensionality in high-dimensional data can be vital in suggesting such 
observables from simulations. Principal Component Analysis is an established linear 
method for the detection of appropriate lower-order observables from simulation data; 
numerically estimated eigenmodes of the problem may also provide good observables 
(see the discussion in jSH] about estimating gaps between eigenvalues, and using them 
to decide whether we should include more observables as independent variables). There 
are, however, some important developments in this area: the recent use of harmonic 
analysis (geometric diffusion) on graphs constructed from high-dimensional data shows 
great promise in detecting good observables (reaction coordinates) for complex, high 
dimensional systems (201111 Hill • This "variable-free" approach can be naturally linked 
to equation-free computation (one designs computational experiments both to detect 
the appropriate observables and to do computations with them) we are currently 
working on demonstrating this link for gene regulatory network modeling. 

It is clear that, in certain cases, an equation-free computational approach is expected 
to have advantages over direct simulation. For steep potentials and low noise, for exam- 
ple, the way equation-free computation uses a good observable to bias the simulation 
will sample the effective potential and give a good estimate of the transition rates much 
faster than direct simulation. Also, parametric analysis methods should be able to ex- 
plore parametric transitions faster and more systematically than direct simulation, in 
analogy with the use of bifurcation techniques rather than direct simulation in deter- 
ministic dynamical systems (e.g. by writing augmented algorithms that converge on 
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marginally stable or unstable solutions). The complexity of the computation depends 
crucially on the dimensionality of the unavailable reduced model, and not so crucially on 
the dimensionality of the detailed, full model. We are currently working on the quan- 
tification of these computational benefits; this work is complicated by the fact that - 
lacking explicit formulas from which to obtain derivative information - errors must be 
computed on-line, through a posteriori estimates. 

This brings us to a final, yet vital issue: estimation. Given the noisy nature of the 
data, estimating the numerical quantities of interest lies at the heart of the accuracy (and 
thus the viability) of the computation. For our gene networks, these quantities included 
the effective potential <&(q) and the effective diffusion coefficient D(q). Preliminary in- 
vestigations revealed that the effective diffusion coefficient, D(q), is quite sensitive to the 
time step At and needs to be treated with care. Also, small changes in the average ve- 
locity or effective diffusion coefficient can have relatively large effects on the steady state 
distribution. Even though some computations are "embarrassingly parallel" (one short, 
fine scale realization per processor, running independently) variance reduction becomes 
an important feature (see, e.g. [2E])- Maximum likelihood estimation techniques (e.g. 
[2]) take the place of simple formulas such as ()3.2j) and (J3.3)) : one can envision certain 
hypothesis testing computations (is our model locally well-approximated by a diffusion 
process ?) becoming part of the overall computational scheme. Until these elements, 
and their computational cost, are analyzed and tested, there will be no firm guarantees 
for the computational efficiency of equation- free methods. Yet, even with these caveats, 
as we computationally demonstrated in this paper, we believe that the equation-free 
framework provides a promising new approach to gene regulatory network modeling, al- 
ternative to long-direct simulation. It links directly with powerful and tested traditional 
continuum numerical algorithms (such as numerical integration, fixed point algorithms, 
matrix-free iterative linear algebra) and with system theory techniques like filtering and 
estimation. These techniques are, in some sense "off the shelf and do not need to be 
redeveloped. In our opinion, it is the linking of equation-free techniques with novel data 
reduction/clustering techniques (such as the use of the graph Laplacian to detect good 
reaction coordinates [SU]) that hold the most promise in the computational study of 
complicated stochastic systems in general, and of gene regulatory networks and their 
models in particular. 
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